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1.1 Introduction 

In principle, we should not need the time-dependent extension of density- 
functional theory (TDDFT) for excitations, and in particular not for Molecu- 
lar Dynamics (MD) studies: the theorem by Hohenberg and Kohn [Hohenberg 
1964] teaches us that for any observable that we wish to look at (including 
dynamical properties or observables dependent on excited states) there is a 
corresponding functional of the ground-state density. Yet the unavailability 
of such magic functionals in many cases (the theorem is a non-constructive 
existence result) demands the development and use of the alternative exact 
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reformulation of quantum mechanics provided by TDDFT. This theory de- 
fines a convenient route to electronic excitations and to the dynamics of a 
many-electron system subject to an arbitrary time-dependent perturbation 
(discussed in previous chafers of this book). This is, in fact, the main pur- 
pose of inscribing TDDFT in a MD framework -the inclusion of the effect of 
electronic excited states in the dynamics. However, as we will show in this 
chapter, it may not be the only use of TDDFT in this context. 

The term "Ab Initio Molecular Dynamics" (AIMD) has been exclusively 
identified in the past with the Car-Parrinello (CP) technique [Car 1985]. This 
method combines ground-state DFT with MD, providing an efficient refor- 
mulation of ground-state Born-Oppenheimer MD (gsBOMD) [Marx 2000]. 
However, the "AIMD" words have broader meaning, and should include all 
the possible MD techniques that make use of a first principles approach to 
tackle the many-electron problem. For example, Ehrenfest MD (EMD) can 
also be one AIMD scheme if TDDFT is used to propagate the electronic sub- 
system. This is the most common manner in which TDDFT and MD have 
been combined in the past: as a means to study fast out-of-equilibrium pro- 
cesses, typically intense laser irradiations or ionic collisions [Saalmann 1996, 
Saalmann 1998, Reinhard 1999, Kunert 2001, Castro 2004]. 

Nevertheless, there are other possibilities. In this chapter we review two 
recent proposals: In Section II. 2[ we show how TDDFT can be used to de- 
sign efficient gsBOMD algorithms [Alonso 2008, Andrade 2009] -even if the 
electronic excited states are in this case not relevant. The work described 
in Section 11.31 addresses the problem of mixed quantum-classical systems at 
thermal equilibrium [Alonso 2010]. 



1.2 Fast Ehrenfest IVIolecular Dynamics 

Ehrenfest Molecular Dynamics (EMD) is a model for describing the evolution 
of a mixed quantum-classical system. The equations of motion are given by: 

M^^Rait) = -{m\^RMRit),t)Mt)) . (1-1) 

i^^mt)) ^HiRit),mit)), (1-2) 

where ^(t) is the state of the quantum subsystem (we will assume that this 
is a set of N electrons), and {Ra}^^{; are the position coordinates of A^nuc 
classical particles (a set of N^uc nuclei of masses Ma and charges Za). The 
quantum (or electronic) Hamiltonian operator He{R,t) depends on these 
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classical coordinates, and is usually given by: 
^ 1 1 

t=l i,3<i ' f}<a ' 

- E \R^-rA + E <xt (r. , i) + E ^ext {Rc.t), (1.3) 

Q.z i a 

where v°^^. and w^^t are external potentials acting on the electrons and nu- 
clei, respectively [Echenique 2007]. Atomic units are used throughout the 
document in order to get rid of constant factors such as h or I/A-keq. 

Given this definition, one can show that Eq. (jl.ip can be rewritten as: 

M^—R^{t) = - /dV n(r, t)VR„«o(r, R{t)) , (1.4) 

where 

a /3<a ' P' 

a result which is known as the "electrostatic force theorem" in the quantum 
chemistry literature [Levine 2000], and which is based in the fact that the 
gradient \7 ji^Hc{R{t),t) is a one-body local multiplicative operator (as far 
as the electrons are concerned), i.e., it is a sum of one-electron operators 
whose action amounts to a multiplication in real space [Eschrig 2003, Von 
Barth 2004]. 

Eq. (|1.4p shows that the knowledge of the time-dependent electronic den- 
sity n{r,t) suffices to obtain the nuclear movement. This fact is the basis for 
TDDFT-based Ehrenfest MD (E-TDDFT): instead of solving Eq. ([Tl]), we 
solve the corresponding time-dependent Kohn-Sham system, which provides 
an approximation to n{r,t): 

i^^iPj{r,t) = -^V^iPjir,t)+VKs[n]{r,t)ip,ir,t), j^l,...,N, (1.6) 
being 

VKs[n]ir,t) := E i p 7^^" : + vn[n]{r,t) + v^,[n]{r,t) + vl,,{r,t) , (1.7) 

^ \Ra(t) - r 

and 

N 

n(r,t):==2El^.(r-,i)P, (1-8) 
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where vks [fA {f'l t) is the time-dependent Kohn-Sham potential, and uh [n] (r , t) 
and Vxc[n]{r,t) are the Hartree and exchange-correlation potential, respec- 
tively. For simplicity, we assume an even number of electrons in a spin- 
compensated configuration. 

The equations of motion for E-TDDFT ([131 [H] and HHl) can be derived 
from the following Lagrangian (assuming an adiabatic approximation for the 
exchange and correlation potential, as it is commonly done in practical im- 
plementations of TDDFT): 



for /X = 1 (the reason for including this parameter fi will become clear in 
what follows). We use a dot to denote time-derivatives. 

The term i?KS is the Kohn-Sham ground-state energy functional: 



Note that, when the time-dependent orbitals are introduced into this ex- 
pression (as it is done in E-TDDFT), it becomes a functional of the Kohn- 
Sham orbitals at each time, and not a functional of the ground state density. 
Also, from here on, we assume that there are no external potentials Vp^t and 
I'gxt) since they do not add anything to the following discussion. 

It is worth remarking now that EMD differs from gsBOMD, and it is in- 
structive to see in which way. We do so in the initial formulation of Eqs. (jl.ip 
and (|1.2p using the A'^-electron wavefunction for simplicity, i.e., we forget for 
a moment the TDDFT formalism. 

To illustrate the main concepts we start by projecting the EMD equations 
into the adiabatic basis, formed at each nuclear configuration by the set of 
eigenfunctions of the electronic Hamiltonian: 




(1.9) 



a 





(1.10) 



H,{R)\^„,{R)) = E„,{RMn{R)) , 



(1.11) 




(1.12) 



m 
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The result is: 

il 



(1.14) 



dt 

where the "non-adiabatic couplings" are defined as: 

dr(i«) := {^m{R)\'7RMR)) ■ (1-15) 

If these are negligible, and we assume that the electronic system starts 
from the ground state (0^(0) = 6mo), EMD reduces to gsBOMD: 

d2 

Mo,—Ra{t) = VnMRit)) : (1-16) 

Cmit) =SmO- (1-17) 

Now, in order to integrate the gsBOMD equations, one can make use of 
ground-state DFT, since the only necessary ingredient is the ground-state 
energy EQ{R{t)). One could thus prccomputc this hypcr-surfacc, in order to 
propagate the nuclear dynamics a posteriori, or else only compute the energies 
at the -R points visited by the dynamics (a procedure normally known as "on- 
thc-fly"). However, Car and Parrinello [Car 1985] proposed an alternative, 
based on the following Lagrangian: 



LJP[VP, if, R, it] := $^(<^j#,> + J2 \^ocRI 

3 a 

-^Ksb, R] + l^^ij - ■ (1.18) 



Note the presence of a fictitious mass A, and of a set of Lagrange multipli- 
ers Aij associated to the constraints that keep the KS orbitals orthonormal 
along the evolution. The Car-Parrinello (CP) equations that stem from this 
Lagrangian are: 

M^^Ro,{t) = -VR„E^sW),R{t))\ , (1.19) 
\(pj{r,t) = -^S7'^ipj{r,t) + VKs[n]{r,t)ipj{r,t) + ^AjkMr,t) , (1.20) 

k 

{Mt)\Vjit)) = 5ij . (1.21) 
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The first of these three sets of equation ensures that CP molecular dynam- 
ics (CPMD) is (approximately) equivalent to gsBOMD if the KS orbitals stay 
close to the ground-state ones; the second equation is an auxiliary, fictitious 
electronic propagation that enforces this proximity to the ground state for a 
certain range of values of the "mass" A; whereas the last equation demands 
the constant orthonormality of the electronic orbitals. Another role of the 
fictitious mass lambda is to accelerate the fake electronic dynamics, and as 
a consequence to improve the numerical efficiency. This efficiency (in addi- 
tion to the success of DFT in the calculation of total energies with chemical 
accuracy) has made of CPMD the method of choice for performing ab initio 
gsBOMD during the last decades. 

When attempting simulations of very large systems, the calculations must 
be done using the massive parallel architectures presently available, therefore 
one must ensure a good scalability of the computational algorithms with 
respect to the number of processors and the size of the systems (ie. number 
of atoms). The CPMD technique at a given point has to face the problem 
posed by the need of orthonormalization, as required by Eq. (jl.2ip . This is a 
very non-local process (regardless of the algorithm used) , and therefore very 
difficult to parallelize efficiently. Linear-scaling methods and other approaches 
have been proposed recently [Kiihne 2007] to improve the speed of the CP 
technique. 

One possibility to circumvent the orthonormalization issue is to do E- 
TDDFT (which automatically conserves the orthonormality) instead of CPMD, 
for those cases in which the coupling to higher electronic excited states is 
weak, and therefore E-TDDFT is almost equivalent to gsBOMD. This fact 
was first realized by Theilhaber [Theilhaber 1992]. Unfortunately, the re- 
quired time step for E-TDDFT is very small (two to three orders of magnitude 
smaller than the CPMD time-step), which makes it very inefficient compu- 
tationally. The reason is that the simulation must follow the real electronic 
motion, which is very fast (in contrast to the fictitious electronic motion used 
in CPMD). In [Alonso 2008] and [Andrade 2009], however, it was shown how 
the time-step can be increased by modifying the parameter in the defini- 
tion of the Lagrangian function given in Eq. (jl.9p . which for normal Ehrenfest 
dynamics should he ^ — 1. 

For any /i, the equations of motion derived from this Lagrangian function 

are: 

ivVj(r-,i)+WKsW(r,t)(^j(r,t), (1.22) 

VR^EKsW),R{t)]- (1-23) 

The only difference with respect to E-TDDFT is the appearance of the 
jjL parameter multiplying the time-derivative of the time-dependent KS equa- 
tions. The most relevant features of this dynamics are: 



d^ 

M^—R^{t) = - 
at^ 
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1. The orthogonality of the time-dependent KS orbitals is automatically 
preserved along the evolution, so that there is no need to perform any 
orthonormalization procedure. 

2. The "exact" total energy of the system, defined as 

^= ^E^-^"^' +^Ks[^,H], (1.24) 

is also preserved along the evolution. Note that it is independent of /x 
and it coincides with the same exact energy that is preserved along the 
gsBOMD evolution. In contrast, the preserved energy in CPMD is given 
by: 

EcP = E + ]^\Y,{^iW3) ^ (1.25) 

3 

where E, given by Eq. (I1.24p . is now time-dependent. It can be seen how 
the new constant of motion £'cp actually depends on A, which is the 
fictitious electronic mass introduced in the CP formulation. 

3. If we consider Eq. (|1.22p . and write the left hand side as: 

. d d , , 

the resulting equations can be seen as the standard Ehrenfest method 
in terms of a fictions time i^j. This has the effect of scaling the TDDFT 
excitation energies by a l//x factor. So we may open or close the electronic 
gap by using a ^ smaller or larger than 1. Obviously, if /i — 0, then the 
gap becomes infinite, and we retrieve the adiabatic (gsBOMD) regime. 

4. The second important effect of this time re-scaling is a change in the 
required time-step for the numerical propagation; if the time-step for the 
standard E-TDDFT equations is At, then the required time-step for the 
new dynamics is At^^ = fiAt. In other words, the propagation will be /i 
times faster. 

5. Taking into account the two previous points and recalling that the pur- 
pose of this modified Ehrenfest dynamics is to reproduce, albeit approx- 
imately, the gsBOMD results, it becomes clear that there is a tradeoff 
affecting the optimal choice for the value of ^: low values (but still larger 
than one) will give physical accuracy, while large values will produce a 
faster propagation. The optimal value is the maximum value that still 
keeps the system near the adiabatic regime. It is reasonable to expect 
that this value will be given by the ratio between the electronic gap and 
the highest vibrational frequency of the nuclei. For many systems, like 
some molecules or insulators, this ratio is large and we can expect large 
improvements with respect to standard Ehrenfest MD. For other systems, 
like metals, this ratio is small or zero and the new method will not work. 
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Number of atoms Number of processors 

Fig. 1.1. (a) Scheme of the benzene molecule array, (b) Single processor compu- 
tational cost for different system sizes, (inset) Polynomial extrapolation for larger 
systems, (c) Parallel computational cost for different system sizes, (d) Parallel scal- 
ing with respect to the number of processor for a system of 480 atoms. In both 
cases, a mixed states-domain parallelization is used to maximize the performance. 

6. Regarding the scaling with the system size, the modified Ehrenfest dy- 
namics evidently inherits the main advantage of the original one: since 
propagation preserves the orthonormality of the KS orbitals, it needs not 
be imposed and the numerical cost is proportional to NwNc (with Nw 
the number of orbitals and Nc the number of grid points or basis set 
coefficients). For CPMD, a reorthogonalization has to be done each time 
step, so the cost is proportional to N^Nc- From these scaling properties, 
we can predict that for large enough systems the Ehrenfest method will 
be less costly than CP. For smaller systems, however, this gain will not 
compensate for the fact that the time-step, despite being increased by 
the fi factor, will still need to be one or two orders of magnitude smaller 
then the time-step utilized in CPMD. 

Some numerical examples, performed with the octopus code [Marques 
2003, Castro 2006], that give an idea of the performance of this modified 
Ehrenfest dynamics were shown in [Alonso 2008] and [Andrade 2009] . We re- 
produce here one case: the vibrational spectrum of an artificial benzene cristal 
(see Fig. Essentially, the calculations consist of the time-propagation 

of the system, either with the standard CP technique or with the modi- 
fied Ehrenfest dynamics, for an interval of time departing from a Boltzmann 
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distribribution of velocities at a given temperature. Then, the vibrational 
frequencies are obtained from the Fourier transform of the velocity autocor- 
relation function. 

Panels (b) and (c) display the serial and parallel computational cost, 
respectively, of the two methods, defined as the computer time needed to 
propagate one atomic unit of time. In the serial case, it can be seen how for 
the system sizes studied, CP is more efficient; a different scaling can already 
be guessed from the curve; indeed, if these curves are extrapolated (inset), 
one can predict a crossing point where the new Ehrenfest technique starts 
to be advantageous. This is more patent in the parallel case, as can be seen 
in panel (c). Panel (d) displays the different scalability of the two methods: 
for a fixed system size, the system is equally divided in a variable number of 
processors, and the figure displays the different speed-ups obtained. 

The key conclusion is that the lack of the orthonormalization step per- 
mits a new efficient parallelization layer, on top of the usual ones that are 
commonly employed in CPMD (domain decomposition, and Brillouin zone 
K-points): since the propagation step is independent for each orbital, it is 
natural to parallelize the problem by distributing the KS states among pro- 
cessors. Communication is only required once per time-step to calculate quan- 
tities that depend on a sum over states: the time dependent densities and the 
forces over the ions. 



1.3 MD at finite electronic temperature 

The previous section has addressed algorithmic alternatives to the solution 
of the gsBOMD equations ()1.16|) and (|1.17p . These represent the evolution 
of the classical nuclei, interacting all-to-all through the potential Eo{R). The 
resulting dynamics can be used to calculate equilibrium averages at a given 
finite temperature, by assuming ergodicity and computing time averages over 
a number of trajectories, once the system has been appropriately coupled 
to a thermostat. The resulting marginal equilibrium density in the nuclear 
positions space being, in the canonical ensemble: 

with /3 :— l/fcsT, or (3 1/RT if per-mole units are used. 

However, this scheme ignores completely the dynamics of the electrons, 
by assuming that, even at a finite temperature, they are continuously tied 
to their ground state. This assumption is legitimate if the electronic gap 
is large compared to kgT at the temperature of interest. Indeed, in many 
physical, chemical or biological processes the dynamical effects arising from 
the presence of low lying electronic excited states have to be taken into ac- 
count. For instance, in situations where the Hydrogen bond is weak, different 
states come close to each other and non-adiabatic proton transfer transitions 
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become rather likely at normal temperature [May 2004]. In these circum- 
stances, the computation of ensemble averages cannot be based on a model 
that assumes the nuclei moving on the ground-state BO surface. 

In the DFT realm, the inclusion of electronic excited states in the dy- 
namics is very often done by working with partial occupation numbers to 
account for the electronic excitations [Grumbach 1994, Alavi 1994, Alavi 
1995, Marzari 1997], ideally making use of temperature-dependent exchange 
and correlation functionals [Mermin 1965, Prodan 2010, Eschrig 2010]. This 
scheme is however tied to DFT, and is hindered by the difficulty of realisti- 
cally approximating this functional. Other alternative options are Ehrenfest 
dynamics and surface hopping [Tully 1990] (for more on recent progress in 
non-adiabatic electronic dynamics in mixed quantum-classical dynamics, see, 
for example, [Zhu 2005]). 

Recently, Alonso et al. [Alonso 2010] have proposed a new alternative, 
which can make use of the ability of TDDFT to compute electronic excited 
states. In the following, we make a summary of the new technique. 

In order to arrive to a general quantum-classical formalism, and to a 
suitable expression for the quantum-classical equilibrium distribution that is 
considered to be the correct one in the literature, it is preferable in this case 
to follow the partial Wigner transformation route [Wigner 1932], as done 
by Ciccotti, Kapral and Nielsen [Kapral 1999, Nielsen 2001]. Let us assume 
a quantum system of two particles of masses m and M (M > m) living 
both in one dimension, whose canonical position and momentum operators 
are and {X,P), respectively. The generalization to more particles and 

higher dimension is straightforward. Given an operator A, its partial Wigner 
transform Aw with respect to the large-mass coordinate is defined as: 

Aw{X, P) := (27r?l)-i Jdz e'^^/^(X - z/2\A\X + z/2) . (1.28) 

The operator Aw (X, P) acts on the Hilbert space of the light particle, 
and depends on the two real numbers (X, P). It is possible to reformulate all 
quantum theory in terms of these partial Wigner transforms; in particular, if 
the Hamiltonian for the two particles is given by: 

its transformation is: 

Hw{X,P) = ^ + l^ + V[x,X), (1.30) 

i.e., one just has to substitute the quantum operators of the heavy particle 
by the real numbers (AT, P) . 

If the state of the system is described by the density matrix p{t), its 
evolution will be governed by von Neumann's equation, 

d , i r_^_ ,1 , _ , 
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which can be cast into its partial Wigner-transformcd form: 



hA/2i 



Pw - Pw 



w 



where A is the "Poisson bracket operator" , 



A 



d d 
'dPdX 



d d 
dXdP 



(1.32) 



(1.33) 



and the arrows indicate the direction in which each derivative acts. 

Note that up to now, this is an exact reformulation of quantum mechanics 
(no classical or semiclassical limit has been taken). However, this is also a 
convenient departure point to take the classical limit for the heavy particle. 
After an appropriate change of coordinates [Kapral 1999], if we retain only 
the first order terms in r] := (m/My^^, Eq. (|1.32p is transformed into: 



d . t 



Pw 



2 ( {Hw, Pw] — {pwi Hw 



(1.34) 



where {•,•} is the Poisson bracket with respect to the canonical conjugate 
coordinates {X,P), 



dHw dpw dHw dpw 
dX ~dP dP~~dx' 



(1.35) 



and both pw and Hw are functions of {X, P). 

The equilibrium density matrix in the partial Wigner representation at 
the classical limit for the heavy particle, denoted by p^ should be stationary 
with respect to the evolution at first order in r] := {m/MY^'^ in Eq. (|1.34|) . 
If we use this property and expand the equilibrium density matrix in powers 
of 7]-. 



if^{X, P)^^ i^-ff^ ^"^ (X, P) , (1.36) 

ra=0 

it can then be proved [Nielsen 2001] that the zero-th order term is given by: 
p;^(°)(X,P) = le^^^«-(^^^), (1.37) 

with 



Z := Tra 



dXdPe-^^-(^^^) 



(1.38) 



the symbol Trq meaning trace over the quantum degrees of freedom. 

Note that (I1.37P corresponds, at fixed classical variables (X, P), to the 
equilibrium density matrix for the electronic states. However, it is only an ap- 
proximation to the true quantum-classical equilibrium density matrix, since 
it is not a stationary solution to the quantum-classical Liouvillian given 
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in Eq. (I1.34p . This distribution is often regarded, however, as the correct 
equilibrium distribution of the canonical ensemble for a mixed quantum- 
classical system [Mauri 1993, Parandekar 2005, Parandekar 2006, Schmidt 
2008, Bastida 2007], and the average of observables is computed as: 

{dix,p,X,P)) ^Tr^J dXdP6(x,p,X,P)p°^(°)(X,P). (1.39) 

As mentioned, the careful analysis described in [Kapral 1999, Nielsen 
2001], shows that this is a first order approximation in the square root of the 
quantum-classical mass ratio rj := {m/My/^, and therefore an acceptable 
approximation if this ratio is small. 

In the remaining part of this chapter, and following [Alonso 2010], we 
will write a system of dynamic equations for the classical particles such that 
the equilibrium distribution in the space of classical variables is in fact given 
by Eq. (I1.37p . This is also a goal of surface hopping methods [TuUy 1990], 
although it is not fully achieved since these methods do not exactly yield 
this distribution [Schmidt 2008]. We will do this by deriving a temperature- 
dependent effective potential for the classical variables, which differs from the 
ground-state potential energy surface (PES) used in gsBOMD. It is straight- 
forward, however, to write an equation that gives the expression for the effec- 
tive potential in terms of this PES together with the BO PESs corresponding 
to the excited states of the electronic Hamiltonian. Despite this property, it 
is worth remarking that the approach described here is based on the assump- 
tion that the full system of electrons and nuclei is in thermal equilibrium at 
a given temperature, and not on the assumption that electrons immediately 
follow the nuclear motion (i.e., the adiabatic approximation), which is at the 
core of the BO scheme. 

Let us assume that we are only interested in the average of observables 
that depend explicitly only on the degrees of freedom of the heavy, classical 
particle, A — A{X, P). It is a matter of algebra (using Eqs. (|1.37|) and ()1.39p ) 
to prove that this average can be written as: 




dXdPA{X,P)c-'^"'="'^^-^-^f''> , (1.40) 



where we have introduced an effective Hamiltonian H^s, defined as: 

H,n{X,P-p) := -llnTrqC-^^-^^'^) . (1.41) 

The partition function Z can also be written in terms of the effective 
Hamiltonian: 

Z - y dXdPe-'^^-(^'^^'" , (1.42) 

Hence, the quantum subsystem has been "integrated out" , and does not 
appear explicitly in the equations any more (of course, it has not disappeared. 
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being hidden in the definition of the effective Hamiltonian) . In this way, 
the more complicated quantum-classical calculations have been reduced to a 
simpler classical dynamics with an appropriate effective Hamiltonian, which 
produces the same equilibrium averages of classical observables [Eq. (|1.40|) ] 
as the one we would obtain using Eq. (|1.37p in ()1.39|) . and hence incorporates 
the quantum back-reaction on the evolution of the classical variables, at least 
at the level of equilibrium properties. 

In the case of a molecular system, the total (partially Wigner transformed) 
Hamiltonian reads: 



where R denotes collectively all nuclear coordinates, P all nuclear momenta, 
T„(P) is the total nuclear kinetic energy, and Hc{R) is the electronic Hamil- 
tonian in Eq. (|1.3|) . that includes the electronic kinetic term and all the 
interactions. The effective Hamiltonian, defined in Eq. (|1.41|) in general, is in 
this case of a molecular system given by: 



ifeff(H,P;/3) :=rn(P) - -InTrqC^'^^^t-^) =: T^iP) + V,s{R; /3) , (1.44) 



where the last equality is a definition for the effective potential Vcs{R] /3). 

Now, making use of the adiabatic basis, defined in Eq. (|l.lll) as the set of 
all eigenvectors of electronic Hamiltonian Hc{R), we can rewrite Vcs{R;/3) 
as: 



where Eno{R) :— £'„(!?)— -Eo(-R)- It is for the computation of these excitation 
energies that TDDFT can be employed. The proposed dynamics would be, 
therefore, based on TDDFT. Of course, any other many-electron technique 
can also be used. 

This equation permits to see explicitly how the ground state energy Eq 
differs from VcS , and in consequence how a MD based on VcS is going to differ 
from a gsBOMD. In particular, notice that Vcs{R; (3) < Eo{R), and compare 
the marginal probability density in the gsBOMD case in Eq. (|1.27p to the 
one produced using the new dynamics: 



H{R,P)^Tn{P) + H,{R), 



(1.43) 



V,s{R; (3) ^ Eo{R) - -In 1 + ^ e'^^-^^) 



(1.45) 



n>0 



Pes{R) 




(1.46) 



/ dR' (l + ^e 

V Tl>0 



■l3E„o{R') g-/3Bo(H') 



Finally, note that to the extent that nuclei do not have quantum behavior 
near conical intersections or spin crossings, nothing prevents us to use this 
equation also in these cases. 
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The definition of tlie classical, effective Hamiltonian for the nuclear coordi- 
nates in Eq. (II. 44^ allows us now to use any of the well-established techniques 
available for computing canonical equilibrium averages in a classical system. 
Of course, since Hes in Eq. (|1.44p depends on T, any Monte Carlo or dynam- 
ical method must be performed at the same T that -ffeff was computed in 
order to produce consistent results, given in this case by the convenient ex- 
pression (ll.40p . For example, we could use (classical) Monte Carlo methods, 
or, if we want to perform MD simulations, we could propagate the stochastic 
Langevin dynamics associated to the Hamiltonian ()1.44|) : 



MjRj{t) = ~VjVMR{t);f3) - M.nRj{t) + MjS{t) , (1.47) 



where S is a vector of stochastic fluctuations, obeying {Si{t)) = and 
{Si{ti)Sj{t2)) = 2jkBT6ij6{ti — 12) which relates the dissipation strength 7 
and the temperature T to the fluctuations (fluctuation-dissipation theorem). 

Indeed, it is well-known that this Langevin dynamics is equivalent to the 
Fokker-Planck equation for the probability density W{R^ P) in the classical 
phase space [Van Kampen 2007]: 



Any solution to Eq. ()1.48|) approaches at infinite time a distribution 
Wcq(i?, P) such that 9iW^oq(-R, -P) = 0. This stationary solution is unique 
and equal to the Gibbs distribution, Woq(i?,P) = Z"^ q-I3H^!!{R,P;0) jvan 
Kampen 2007]. Thus, the long-time solutions of Eq. (|1.48p . and hence those 
of Eq. (|1.47|) reproduce the canonical averages in Eq. (|1.40l) . This prop- 
erty, which is also satisfied by other dynamics like the one proposed by 
Nose [Nose 1984, Nose 1991] if the Hcs in Eq. (|1.44|) is used, comes out in a 
very natural way from the present formalism while it is yet unclear of other 
ab initio MD candidates for going beyond gsBOMD [Mauri 1993, Parandekar 
2005, Schmidt 2008, Bastida 2007]. 

When would this new MD scheme be useful? The approach introduced 
in this section is particularly suited to the case of conical intersection or 
spin-crossing [Yarkony 1996], since it does not assume that the electrons or 
quantum variables immediately follow the nuclear motion, in contrast to any 
adiabatic approach. Another interesting application pertains the debated is- 
sue of quantum effects in proton transfer [Iyengar 2008]. It is a matter of 
current debate to what extent protons behave "quantum-like" in biomolec- 
ular systems (e.g. is there any trace of superposition, tunneling or entangle- 
ment in their behavior?). Recently, McKenzie and coworkers [Bothma 2010] 
have carefully examined the issue, and concluded that "tunneling well below 
the barrier only occurs for temperatures less than a temperature Tq which 
is determined by the curvature of the PES at the top of the barrier." In 



dW{R,P;t) 

dt 



= {H,s{R,P;f3),W{R,P;t)} 

+ 7 51 (Pj + MkBTdp,)W{R, P- 1) . 



(1.48) 



,7 
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consequence, the correct determination of this curvature is of paramount im- 
portance. 

The curvature predicted by the temperature-dependent effective potential 
introduced here is smaller than the one corresponding to the ground state 
PES, in the cases in which the quantum excited surfaces approach, at the 
barrier top, the ground state one. Therefore, Tq would be smaller than that 
corresponding to the ground state PES (see Eq. (8) in [Bothma 2010]), and 
hence the conclusion in this reference "that quantum tunneling docs not play 
a significant role in hydrogen transfer in enzymes" is reinforced by the results 
of the new dynamics. 
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